% set_params
% Jones, Philippon, Venkateswaran (2021)
%

%% Compute Steady-State

in = 0 ; % Initial infected fraction
d  = 0 ; % Initial dead fraction 
s  = 1 ; % Initial susceptible fraction
r  = 0 ; % Initial recovered fraction

delta = deltascale*(deltabar + exp(phi * (kappa * in))-1) ;

m1 = 0 ;
M1 = 0 ;
l1 = 1 ;
L1 = (1-d-kappa*in)*l1 ;

mbar1 = 0 ;

chi1t = chibar1*(1-Deltachi*(1-exp(-rhom*mbar1))) ;

c1 = (1-d-kappa*in)*(l1 - chi1t/2*m1^2) ;
C1 = (N-d)*c1 ;

e       = e0 + (1-d)*ec1*c1*C1 + (1-d-kappa*in)*(1-m1)*el1*l1*L1*(1-M1) ;

gamma    = R0*rho/e ; % Transmission rate for given R0

lambdae = ((1/c1)-l1^eta)/(el1+ec1) ;
lambda  = l1^eta + lambdae*el1 ;
lambdad = (1/((1-(1/beta)))) * ( l1^(1+eta) / (1+eta) - ud - log(c1) - lambda*(l1 - chi1t/2 * m1^2 - c1) + lambdae*(ec1*c1*C1 + (1-m1)*el1*l1*L1*(1-M1))) ;
lambdai = (1/((1-rho-delta*kappa-1/beta))) * ( kappa * l1^(1+eta) / (1+eta) - us*kappa - kappa*lambda*(l1 - chi1t/2 * m1^2) + kappa*lambdae*(1-m1)*el1*l1*L1*(1-M1)  - delta*kappa*lambdad ) ;
lambdas = 0 ;

Vmbar = (1/(1-beta))*(1/2)*lambda*(1-d-kappa*in)*chibar1*Deltachi*rhom*exp(-rhom*mbar1)*m1 ;

Vi = -(1/beta)*lambdai ;
Vs = -(1/beta)*lambdas ;
Vd = -(1/beta)*lambdad ;

f  = 0 ;
Vf = 0 ;

a  = 0 ;
Va = 0 ;
abar = amean ;

params = [...
    beta	
    ec1		
    el1		
    chibar1	
    N		
    gamma	
    deltabar
    phi
    kappa	
    rho		
    e0		
    eta		
    rhom	
    Deltachi
    us		
    ud		
    psim
    deltascale
    rhof
    omega1
    omega2
    rhoa
    varphi1
    varphi2
    amean] ;

xinitial = [ ...
    c1		
    l1		
    m1		
    chi1t	
    mbar1	
    Vmbar	
    lambda	
    lambdae	
    lambdas	
    lambdad	
    lambdai	
    Vi		
    Vd		
    Vs		
    e		
    L1		
    C1		
    M1		
    in		
    d		
    s 		 
    delta
    f
    Vf 
    a
    Va
    abar];

% Run SIR model for initial guess to final system
[s_fin, i_fin, r_fin, d_fin] = sir(rho, gamma, delta, kappa, 0, 50, 999, 1) ;

in = i_fin(end)/1000 ;
d = d_fin(end)/1000 ;
s = s_fin(end)/1000 ;
r = r_fin(end)/1000 ;

delta = deltascale*(deltabar + exp(phi * (kappa * in))-1) ;

m1 = 0 ;
M1 = 0 ;
l1 = 1 ;
L1 = (1-d-kappa*in)*l1 ;

mbar1 = 0 ;

chi1t = chibar1*(1-Deltachi*(1-exp(-rhom*mbar1))) ;

c1 = (1-d-kappa*in)*(l1 - chi1t/2*m1^2) ;
C1 = (N-d)*c1 ;

rr = 1/beta - 1 ;
P1 = 1 ;

e       = e0 + (1-d)*ec1*c1*C1 + (1-d-kappa*in)*(1-m1)*el1*l1*L1*(1-M1) ;
lambdae = ((1/c1)-l1^eta)/(el1+ec1) ;
lambda  = l1^eta + lambdae*el1 ;
lambdad = (1/((1-(1/beta)))) * ( l1^(1+eta) / (1+eta) - ud - log(c1) - lambda*(l1 - chi1t/2 * m1^2 - c1) + lambdae*(ec1*c1*C1 + (1-m1)*el1*l1*L1*(1-M1))) ;
lambdai = (1/((1-rho-delta*kappa-1/beta))) * ( kappa * l1^(1+eta) / (1+eta) - us*kappa - kappa*lambda*(l1 - chi1t/2 * m1^2) + kappa*lambdae*(1-m1)*el1*l1*L1*(1-M1)  - delta*kappa*lambdad ) ;
lambdas = 0 ;

Vmbar = (1/(1-beta))*(1/2)*lambda*(1-d-kappa*in)*chibar1*Deltachi*rhom*exp(-rhom*mbar1)*m1 ;

Vi = -(1/beta)*lambdai ;
Vs = -(1/beta)*lambdas ;
Vd = -(1/beta)*lambdad ;

f  = 0 ;
Vf = 0 ;

a  = 0 ;
Va = 0 ;
abar = amean ;

xfinal = [ ...
    c1		
    l1		
    m1		
    chi1t	
    mbar1	
    Vmbar	
    lambda	
    lambdae	
    lambdas	
    lambdad	
    lambdai	
    Vi		
    Vd		
    Vs		
    e		
    L1		
    C1		
    M1		
    in		
    d		
    s 		 
    delta 
    f
    Vf
    a
    Va
    abar];

save input params xinitial xfinal ee_T_i ee_T_vac ee_T_e0 ee_T_d